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ABSTRACT 

Context. Radiative transfer has a strong impact on the collapse and the fragmentation of prestellar dense cores. 
Aims. We present the radiation-hydrodynamics (RHD) solver we designed for the RAMSES code. The method is designed for astro- 
physical purposes, and in particular for protostellar collapse. 

Methods. We present the solver, using the co-moving frame to evaluate the radiative quantities. We use the popular flux-limited 
difi'usion approximation under the grey approximation (one group of photons). The solver is based on the second-order Godunov 
scheme of RAMSES for its hyperbolic part and on an implicit scheme for the radiation diffusion and the coupling between radiation 
and matter. 

Results. We report in detail our methodology to integrate the RHD solver into RAMSES. We successfully test the method in several 
conventional tests. For validation in 3D, we perform calculations of the collapse of an isolated 1 Mq prestellar dense core without 
rotation. We successfully compare the results with previous studies that used different models for radiation and hydrodynamics. 
Conclusions. We have developed a full radiation-hydrodynamics solver in the RAMSES code that handles adaptive mesh refinement 
grids. The method is a combination of an explicit scheme and an implicit scheme accurate to the second-order in space. Our method 
is well suited for star-formation purposes. Results of multidimensional dense-core-collapse calculations with rotation are presented in 
a companion paper. 
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1. Introduction 

Within recent years, star-formation calculations have un- 
dergone a rapid increase in the variety of the physi- 
cal models that are included. The coupling between ra- 
diative transfer and hydrodynamics has been widely stud- 
ied for many year and considering different regime s and 
frames (e.g. M ihalas & Mihalas 1984; Lgw rie et al. 1200 it 
iMihalas & Aueij bOOlFTKrumholz et aimOOTbl) . Radiation hy- 
drodynamics (RHD) methods have been developed in grid- 
based codes (IStone & NormanI [T992t iHaves & NormanI l2003t 



Krumholz et al.' '2007a'; K uiper et alj |201(]|; ISekora & Stond 
2010; Tomida et al . 2010) and also in smoothed particles hy- 
drodynamics (SPH) codes dBoss et al.l2000l:IWhitehouse & Bate! 



drody n 
l2006t [5 



JStamatellos et al.l2007 ). Most of these studies us e the pop- 
ular flux-limited diffusion approxima tion (FLD, e.g. iMinerbol 
I1978H L evermore & Pomraning|ll98TI) approximation to model 
the radiation transport. 
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In star-formation calculations, the easiest method to take 
into account radiative transfer is to use a barotropic ap- 
proximation, which reproduces approximately the thermal be- 
haviour of the gas during the collapse. However, more accu- 
rate RHD calculations show that a barotropic equation of state 
(EOS) ca nnot account for r ealistic cooling and h eating of the 
gas (e.g. iBoss et aL 1 120001: lAttwood et al.l l2009l Commer9on 
et al. in prep, hereafter Paper II). Recently, using radiation - 
magnetohydrodynamics calculations, ICommer9on et al.l (1201 Ol) 
have shown that the barotropic approximation cannot properly 
account for the combined effects of magnetic field and radia- 
tive transfer in the first collapse and in the first core formation. 
On larger scales, radiative transfer has been found to greatly re- 
duce the fragmentation because of the r adiative feedback owing 
to acc retion and protostellar evolution (lBatell20d9t lOffner et al.l 
I2OO9I) . 

In this study, we present a new RHD solver based on the 
FLD approximation, which we i ntegrate in the adaptive mesh 
refinement (AMR) code RAMSES (lTevssieijl2002h . The solver is 
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consistently integrated in the second-order predictor-corrector 
Godunov scheme of RAMSES, which we modify to account for 
the radiative pressure. But we add an implicit solver to handle 
the radiation diffusion and the coupling between matter and ra- 
diation, which involves physical processes on timescales much 
shorter than the hydrodynamical one. The FLD is easier to im- 
plement in an AMR code than more sophisticated methods that 
would require an additional equation on the first moment of 
the ra diative transfer equation (e.g. Ml model, iGonzalez et al.l 
l2007h . The extension to (ideal) MHD flows presents n o partic- 
ular d ifficulties and has already bee n used ([Commercon et al.l 
based on the so lver presented in lFromang et al.l (l2006l and 



Tevssieretal](l2006h . 



The paper is organized as follows: in Sect. 2 we recall the 
RHD equations in the comoving frame we use and we briefly 
present the FLD approximation. The RHD solver for the RAMSES 
code is presented in Sect. 3. In Sect. 4 the method is tested for 
well-known test cases. As a final test, RHD dense core collapse 
calculations with a very high resolution are performed. Section 
5 summarizes our work and the main results andr present our 
assessment of the method's potential. 

2. Radiation liydrodynamics in tlie flux-limited 
diffusion approximation 

2.1. Radiation hydrodynamics in the comoving frame 

We consider the equations governing the evolution of an 
inviscid, radiating fluid, where radiative quantities are esti- 
mated in the comoving frame and are frequency-integrated 
(Mihalas & Mihalas 1981) 



dip + 
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where p is the material density, u is the velocity, P the thermal 
pressure, ctr is the Rosseland mean opacity, F, is the radiative 
flux, E the fluid total energy E - pe + l/2ptr (e is the inter- 
nal specific energy), crp is the Planck opacity, B - B(T) is the 
Planck function, E^ is the radiative energy and is the radia- 
tion pressure. We see that the term crgFi/c acts as a radiative 
force on the material. The material energy lost by emission is 
transferred into radiation, and radiative energy lost by material 
absorption is added to the material. To close this system, we need 
two closure relations: one for the gas and one for the radiation. 
In this work, we only consider an ideal gas closure relation for 
the material: P = {y - l)e - pksT/fimu where y is the specific 
heats ratio, yu is the mean molecular weight, and e = pCyT is 
the gas internal energy. For the radiation, we use the FLD ap- 
proximation to close the systern of m oment equations jMinerbol 
I1978H L evermore &Pomraningll98ll) . In this work, we consider 
only the simplified case of a grey material, where all frequency- 
dependent quantities are integrated over frequency. We cannot 
use a frequency-dependent model for our purpose because of 
CPU limitation. 

In comparison with the laboratory frame formulation. ICastoH 
(Il972h demonstrated that in the comoving frame, an additional 
advective flux of the radiation enthalpy is not taken into account. 
In the dynamic diffusion regime, where the optical depth t >> 1 
and {v/c)t » 1, this radiative flux can dominate the diff'usion 
flux, emissio n or absorption. For an al ternative mixed frame for- 
mulation, see iKrumholz etaP (l2007bh . In the low-mass star do- 



main, the main focus of this work, we do not expect to encounter 
dynamic difiiision situations. 

2.1 .1 . The flux-limited diffusion approximation 

As mentioned earlier, we need a closure relation to solve the 
moment equations coupled to the hydrodynamics (closed by the 
perfect gas relation), and such a relation is of prime importance. 
Many possible choices for the closure relation exist. Among 
these models, the FLD approximation is one of t he simplest ones 
and uses moment models of ra diation transport (lMinerbolll978l: 
ILevermore &Pomraningil981b . 

Under the flux-limited diffusion approximation, the radiative 
flux is expressed directly as a function of the radiative energy 
and is proportional and colinear to the radiative energy gradient 
(Fick's law). Under the grey approximation, we have 



cA 

Fr = V£„ 



(2) 



where A = A(R) is the flux limiter, and R - |V£'rl/(o"R£' r). In this 
study, we retain the flux limiter that has been derived bv lMinerbol 
(Il978h . assuming the intensity as a piecewise linear function of 
solid angle 



A = 



2/(3 + V9 + 12/?2) if 0<R< 3/2, 
{1+R + Vl + 2R)-^ if 3/2 <R< 00. 



(3) 



The flux limiter has the property that A 1/3 in optically 
thick regions and A l/R in optically thin regions. We re- 
cover the proper value for diffusion in the optically thick regime, 
F = -c/(3crR)V£'r, and the flux is limited to ciSr in the optically 
thin regime. Under the FLD approximation, the radiative trans- 
fer equation is then replaced by a unique diffusion equation on 
the radiative energy 



^_V.(— V£,| = c^p(47rB- 
<9f \o-R 



cE,). 



(4) 



3. A multidimensional radiation hydrodynamics 
solver for RAMSES 

3.1. The AMR RAMSES code 



We use the RAMSES code (lTevssieijl2002h . which integrates 
the equations of ide al mag netohydrodynamics (iFrom ang et al] 
I2OO6I: iTevssier et al.l [2006^ using a second-order Godunov fi- 
nite volume scheme. The MHD equations are integrated using 
a MUSCL pred ictor-corrector scheme, originally presented in 
Ivan Leeil(IT97l . Fluxes at the cell interface are estimated with 
an approximated Riemann solver (Lax-Friedrich, HLL, HLLD, 
etc.). For its AMR grid, RAMSES is based on a 'tree-based" 
AMR structure, the refinement is made on a cell-by-cell basis. 
Various refinements can be used (fluid vaiiable gradients, insta- 
bility wavelength, etc.). 

The AMR code RAMSES has often been used for 
star-formation 



purpo ses 



( Hennebelle & Fromang 



2008 



2009 



Hennebelle & Teyssier" '2008'; 'Hennebelle & Ciardil 
Commercon et al. 2010). Commercon et al. (2008) have thor- 
oughly and successfully compared its results with standard 
SPH, showing a good agreement between the methods. 
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3.2. The Eulerian approach and properties of conservation 
laws 

In Eulerian hydrodynamics, the mesh is fixed and gas 
density, velocity, and internal energy are primary variables. 
Eulerian methods fall i nto two groups: finite difiference methods 
(e.g. the ZEUS code, IStone & Normanlll992l: iTurner & Stoiiel 
200 lb and fin ite volume methods (e.g. the RAMSES code, 
Tevssiedl2002l) . In the first group, flow variables are conceived 
as being samples at certain points in space and time. Partial 
derivatives are then computed from these sampled values and 
follow Euler equations. In the finite volume approach, flow 
variables correspond to average values over a finite volume - 
the cell - and obey the conservation laws in the integral form. 
Their evolution is determined through Godunov methods by 
calculating the flux of every conserved quantity across each cell 
interface. 

For an inviscid, compressible flow, the Euler equations in 
their conservative form read 



d,p + V[pu] = 
d,pu + V[pu®u + PI] = 
d,E + V[u(£ + P)] = 



(5) 



This system can be written in the general hyperbolic conserva- 
tive form 



as follows: A - 1/3 + (A - 1/3). We thus distinguish a diffii- 
sive part (Eddington approximation. Pi = l/3ErT) and a cor- 
rection part. The computation of predicted states and fluxes in 
the MUSCL scheme is made under the Eddington approxima- 
tion, which is then corrected in an additional corrective step. The 
RHD equations can be rewritten as 



d,p + V [pu] 

d,pu + W[pu0u + {P+ l/3E,)I] 
d,ET+W['a(Ej + P+l/3E,)] 



d,E, + V [uf:,] 





-pVcD- l/SFEr 
-pu • VO 

-(^- l/3)(uVf:,. + £,V : u) 
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-PrV : u -I- Kppcia^T'^ - E,) 

+v ■(■^^eS 



(8) 

The new system c?,U + VF(I[J) - S (U) is composed of the mod- 
ified hyperbolic left hand side (LHS) and the right hand side 
(RHS) source, corrective and coupling terms S (U) = 5exp+5imp. 
The hyperbolic system, as well as the source and corrective 
terms S exp, are integrated in time with an explicit scheme. The 
modified RHD hyperbolic system reads 



U = 



p 




pu 


pu 


, F(U) = 


pu®u + (P+ l/SEr)! 




u{Et + P+1/3E,) 


■Er. 




UEy 



(9) 



+ v.. 

dt 



<{V) = 0, 



(6) 



where the vector U = (p,pu, E) contains conservative variables, 
and the flux vector F(U) = (pu,pu ® u + PI, u (E + P)) is a linear 
function of U. 

In this paper, we use the second-order Godunov method, but 
applied to the modified Euler equation system under the ELD 
approximation. 

3.3. The conservative radiation hydrodynamics scheme 

Let us rewrite the grey RHD equations under the ELD ap- 
proximation within the comoving frame, taking into account the 
gravity terms 



<9,p + V|jou] = 

d,pu -H V [pu ® u + PI] = -pVO - AVEr 
d,Ej -H V [u (Et + P)] = -pu ■ VcD - P^V : u - AuVE, 
+V ■ ( ^V£,-) 

\PKR V 

d,E, + V [uE,] = -P,V : u + V ■ {^"^E,) 

+Kppc(aTiT'^ - E,) 

(7) 

Note that we rewrite the opacity crj as Kip. The dimension of 
is cm^ g"'. 

The basic idea is to build a solver for a radiative fluid, with 
an additional pressure owing to the radiation field: the radia- 
tive pressure. Following the Euler equations in their conservative 
form, the new conservative quantities are density p, momentum 
pu, total energy Ej of the fluid (gas + photon) per unit volume, 
i.e. pe -I- pu- 12 + E^. Primitive hydrodynamical variables do not 
change for the fluid, but we add a fourth equation for the radia- 
tive energy. 

In order to facilitate these equations in RAMSES and to mini- 
mize the number of changes with the pure hydrodynamical ver- 
sion, we decompose each term where the flux limiter A appears 



This system is used in the predictor-corrector MUSCL temporal 
integration. To predict states, we consider the worst case, where 
the radiative pressure is the greatest (1 /3Ey). For the conserva- 
tive update (corrector step) we consider the LHS of system (jSj. 
The associated eigenvalues corresponding to the three waves are 



Ai = 







4£r 


u - , 




9p 


M 






M + -1 




4E, 


9p 



(10) 



Radiative pressure enlarges the span of solutions, since wave 
speeds are faster Once again, with the Eddington approxima- 
tion, we build the system for the case where the radiative pres- 
sure would be the greatest. Therefore, the waves propagate at a 
speed that is within the wave extrema. 

The next step consists in correcting errors due to the 
Eddington approximation by integrating source terms S ne 



ne — 





-{A-1/3WE, 
-(/I - l/3)(uV£,- + £,-V 
P,-V : u 



u) 



(11) 



We here consider an isotropic radiative pressure tensor P^ = 
AEjI. Other a uthors considered e xtensions to t his closure rela- 
tion using the Levermore (198,4!) FED theory (iTumer & Stonel 
.2001 ; Krumholz et al...2007bi) . 

To ensure the stability of the explicit step, the Courant- 
Friedrich-Levy (CFL) stability condition used to estimate the 
timestep also takes into account the radiative pressure. The up- 
dated CFL condition is simply 



Af < CcFL- 



Ax 



(12) 



u + 
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3.4. The implicit radiative sclieme 

The most demanding step in our time-splitting scheme is to 
deal with the diffusion term V • (^^VE^^ and the coupling term 

Kppc{aj^T'^ - Er), which corresponds to 5 imp. This update has 
to be made with an implicit scheme, because the time scales 
of these processes are much shorter than those of pure hydrody- 
namical processes. Two coupled equations are integrated implic- 
idy 

d,pe ^ -KYpc{a^T^ - E,) 

d,E, - V;5V£r = +Kypc{a^T^ -E,) ' 

which give the implicit scheme on a uniform gricQ 



c„r"+' - c,T" 



At 



Af 

^ cA" 

4p" 



where x is a vector containing radiative energy values. The con- 
jugate gradients (CG) method is one of the most popular non- 
stationary iterative methods for solving large symmetric systems 
of linear equations Ax = b. The CG method can be used if 
the matrix A to be inverted is square, symmetric, and positive- 
definite. The CG is memory-efficient (no matrix storage) and 
runs quickly with sparse matrices. For a N x N matrix, the CG 
converges in less than iterations. Basically, the CG method is 
a steepest-gradient-descent method in which descent directions 
are updated at each iteration. Another advantage is that the CG 
method can be run easily on parallel machines. 

To improve convergence of the CG or even to insure 
convergence if one deals with an ill-conditioned matrix A, we 
use a preconditioning matrix M that approximates A. M is also 
assumed to be symmetric and positive definite. In this work, 
we use a simple diagonal preconditioning matrix, which retains 
only the inverse of A diagonal elements. The convergence of 
the CG algorithm is estimated following two criteria: estimation 



where pe = C,T. The nonUnear term (7""+i)4 makes this 
scheme difficult to invert. Yet it is much easier to solve implicitly 
a linear system. Assuming that changes of temperature are small 
within a time step, we can write 



(14) 

of the norm L?- (criterion 



(T" 

(r+')4 = {T"fli + — 



T") 



'•^"'11 < e) or estimation of the 
norm U" (maximum residual value max{/-'')/max{/"') < e). 
Values of e typically range from 10"^ to 10"^. In appendix we 
present an alternative method to the conjugate gradient, the 
Super- Time Stepping method, which can be used efficiently on 
uniform grids or in some particular cases. 



: 4(r")^r'+' - 3(T")\ 



(15) 

Eventually, with (fT4h). we obtain T"^' as a function of T" and 
E"'^^. Then T"^' can be directly injected in the radiative energy 
equation (fT4b). and E"t^ is finally expressed as a function of 
^"^^'j, E"f\, E". and T". The implicit scheme for the radiative 
energy in a cell of volume V,- in the x-direction becomes 



(£-' -£«,)y, - cAtl^] 5,-.i/2 



+ cAt 



(-) 



r,i /v'-l 



I /-1/2- 



(16) 



//-l/2 Ajc,-_i/2 

= cAt^ .p'l {4a^(T]'fTr' - ia^iT'lf - E^^) 
The gas temperature within a cell is simply given by 



3flR4.cAf (r')' + c,t; + 4,,.cAf£;!ti 
c.. + 4flR4.cAf(r;'f 



(17) 



We compute the Planck and Rosseland opacities and the flux 
limiter with a gas temperature value given before the implicit 
update (with index «) to preserve the linearity of the solver 

3.5. Implicit sclieme integration with the conjugate gradient 
algorithm 

Equation ( fTSI l is solved on a full grid made of cells. It 
results in a system of linear equations, which can be written 
as a linear system of equations 



Ax = b. 



(18) 



' Index n and n+l are used for variables before and after the implicit 
update. Outputs of the explicit hydrodynamics scheme supply variables 
with index n. They do not match the variables at time t" and . 



3.6. Comparison to other schemes 



Other RHD solvers based on the FLD approximation have 
been designed in grid-based codes. Among them, the ZEUS and 
ORION implementations are the most widely used and discussed. 
Compared to ZEU S ( Stone & Norman 1992; Turner &Stoii 
l200lt lHaves"etani2006h . our method is fundamentally different, 
although they also use the comoving frame to estimate the radia- 
tive quantities. ZEUS code is based on a finite difference scheme, 
using artificial viscosity and regular grids. Its non-conservative 
formulation can lead to spurious wave propagation when the res- 
olution is not high enough, or if the radiative pressure dominates 
the characteristic velocity (the classical Burgers equation prob- 
lem). All radiative terms such as the radiation transport and the 
radiative pressure work are integrated implicitly in ZEUS. The 
implicit scheme is based on a Newton-Raphson method, using 
GMRES or LU algorithms for the matrix inversion. 

ORION is a patched-based AM R code, which is les s 
flexible than the tree-b ased AMR (iKrumholz et all l2007ah . 
iKrumholz et alj ( l2007al) implemented the mixed-frame RHD 
equations using a multi-grid, multi-timest ep method to solve the 
impli cit scheme for the radiation module (iHowell & Greenoughl 
l2003h . The hydrodynamics part of ORION uses a second-order 
conservative Godunov scheme, with approximate Riemann 
solvers and very little artificial viscosity to treat shocks and 
discontinuities. Using the same idea as in this study, the dif- 
fusion and matter-radiation coupling terms are integrated im- 
plicitly, while the radiative force and radiativ e pressure work 
are inte grated explicitly. Contrary to our work, IKrumholz et al.l 
(I2007ah do not take into account the radiative pressure in the 
flux estimate at the cell interface for the conservative update, 
which could also lead to an inaccurate wave speed propagation 
in radiation-pressure dominated regions. 
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3. 7. Implicit scheme on an AMR grid 

For studies involving large dynamical ranges, such as star 
formation, it is necessary to extend our implicit scheme to AMR 
grids. The difficulty is to compute the correct fluxes and gradi- 
ents at the interfaces between two cells. We need to carefully 
consider the energy balance on a given volume. Energy balances 
are performed on volumes overlapping two cells, which depends 
on whether the mesh is refined or not. If one considers the face 
of a cell on a level d; three connecting configurations with other 
cells are possible (see Fig.[T]i: 

- Configuration 1: the neighbouring cell is at the same level {: 
cells 1 and 3, 

- Configuration 2; the neighbouring cell is at level £ - I: cells 
1 and i, 

- Configuration 3: 2 neighbouring cells exist at level {+1: cell 
i with cells 1 and 2. 

Last but not least, the difiiision routine is called only once 
per coarse step (no multiple timestepping) and scans the full 
grid, from the finer level to the coarse level. In order to opti- 
mise matrix-vectors products, we choose to avoid dealing with 
configuration 3. Hence, when cells at level { + I are monitored, 
values for cells at level { are updated. Configurations 2 and 3 are 
then performed at the same time. Depending on the configura- 
tion, gradients and flux estimates are diff'erent. In the following, 
we will focus on cell 1, of size Ax x Ax. 



i 


1 


3 








2 


4 



Fig. 1. Example of AMR grid configuration 



3.7.1. Gradient estimate 

Gradients WE,- are estimated between the two neighbouring 
cells centre 



- Configuration 1 : (V£'r)i.3 = 

- Configuration 2 : (VEr)ij - 

3.7.2. Flux estimate 



Er,l - £r,3 

Ax 

Er,\ - Er.i 

3Ax/2 



Let S be the surface of the interface of a cell at level C and 
. the flux across this surface between two cells / and /'. The 
energy rate F xS that is exchanged at this interface is 



- Configuration \\F\^xS' ^ — x S ' 



Ax 



- Configuration 2 : F\ - X S - 



Configuration 3 : f[^^^ x 5^"' = f[^ xS' + fIj^xS*^ 



3 Ax 12 

.t 
/.I 



xS' 



Because access to neighbouring finer cells is not straightfor- 
ward, we see from configuration 3 all the interest of updating 
quantities at level t - \ when scanning grid at level t. 

3.8. Limits of the methods 

The first drawback is the use of the FLD approximation 
that implies isotropy of the radiation field. Anisotropics in the 
transparent regime are not well processe d with the FLD, con- 
trary to more accurate models like Ml dGonzalez et all l2007b 
or VETF (Haves & Norman 2003)). A second limitation comes 
from the grey opacity assumpti on, which could limi t the accre- 
tion on the protostars (see IZinnecker & Yorke 2007[ and refer- 
ences therein). In high-mass star formation, Yorke & Sonnhalterl 
(l2002l) show that using a frequency-dependent radiative transfer 
model enhances the flashlight effect and helps to accrete more 
mass onto the central protostar 

From a technical point of view, our method works only for 
unique time stepping, i.e. all levels evolve with the same time 
step. We do not take advantage of the multiple time stepping 
possibility. As a compromise, we investigate the possibility to 
evolve finer levels with their own time steps and perform a 
diffusion-coupling step every 2, 4 or more finer time steps. As 
a result, we find that performing the diffusion step only every 
2 or 4 fine time steps gives correct results. The frequency of 
the implicit solver calls is left to the user convenience, by use 
of the mutli-time stepping of RAMSES. For instance, for a grid 
of levels ranging form level (mm to ^niax, a unique timestep can 
be used for levels ranging form level to £/. In that case, 
only the levels finer that will use not updated radiative quan- 
tities in the Godunov solver A future development would be to 
use a multigrid solver or pre conditioner for parabolic equations 
dHowell & Greenoughll2003l) . 

Another difficulty comes from the residual norm and scalar 
estimates in the CG algorithm. For a large grid with a large 
number of cells, the dot product can be dominated by round-off 
errors, owing to estimates close to machine precision. This be- 
comes even worse in parallel calculations. The usual MPI func- 
tion MPLSUM fails with a large number of processors, the re- 
sults of any sum becoming a function of number of processors. 
This dramatically affects the number of iterations. We imple- 
mented a new MPI functio n that performs sum mation in double- 
double precision following iHe & Dind (l2000h . using the iKnutl] 
(II997h method. 

Eventually, our method involves only immediate neighbour- 
ing cells, whatever their refinement level. As a consequence, 
our method is only first order at the border between levels. 
This could give rise to a loss of accuracy in diffusion prob- 
lems, because gradient estimates are not second-order accurate 
when neighbo uring cells are at finer levels (see configuration 3, 
lPopinelll2003l) . However, the tests we performed ascertain that 
the method is still roughly second-order accurate. The errors are 
only confined to surfaces much smaller than the total volume. 



4. Radiation hydrodynamics solver tests 

4.1. 1D test: linear diffusion 

We only consider the radiative energy diffusion equation in 
this test, without either hydrodynamics or coupling with the gas. 
The equation to integrate is simply 



'df 



3pKR 



-V£, 



0. 



(19) 
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N=16 




Table 1. CPU time, total number of time steps, and number of 
iterations per time step for various numbers of cells N. 
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Fig. 2. (a): Comparison between numerical solution (squares) 
and analytical solution (red line) at time t=l x 10 for the cal- 
culations with 16 cells, (b): Li norm of the error as a function of 
h = I /Ax. The dotted line shows the evolution of the error as a 
function of and the dashed line the evolution of the error as a 
function of li. 



Consider a box of length L-l. The initial radiative energy cor- 
responds to a delta function, namely it is equal to 1 everywhere 
in the box, except at the centre, where it equals Ei-^/i^x - Eq - 
1 X 10^. To simplify, we choose p/CR = 1 and a constant time step. 
We apply Von Neumann boundary conditions, i.e. zero-gradient. 
The analytical solution in a /^-dimensional problem is given by 



(20) 



where = c/{3pKfd. 

Figure|2ja) shows results at time t - 1x10"'^ for a resolution 
of N - 16 cells. The numerical solution is very close to the 
analytical one, even with this small number of cells. In Fig.|2lb) 
we show the evolution of the Li norm of the relative error as a 
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function of h - I /Ax. The Li norm is estimated as 



Zi \E,j-E,Jxi,t)\Axi 

2f£r,«(X/,f)Ax; 



(21) 



where Ei-j is the numerical value of the radiative energy at posi- 
tion Xi and time f, and E^ aixi, t) the corresponding analytic value. 
The error clearly grows as (dotted line), which indicates that 
our method is second-order accurate. 

In Table[T]we report the CPU time, the total number of iter- 
ations, and the number of time steps for various numerical res- 
olutions. At low resolution, the number of time steps increases 
linearly with the number of cells, as well as the CPU time. The 
number of iterations per time step is constant, i.e. the conver- 
gence of CG does not depend on the dimension of the problem, 
but on the nature of the problem. 

4.2. 1D test: nonlinear diffusion 

In this second test, we consider an initial discontinuity in a 
box with different initial radiative energy states: E^ = 4 on the 
left and E^ - 0.5 on the right. We apply Von Neumann bound- 
ary conditions. We integrate the same equation as in the previous 
test, but with a Rosseland opacity as a nonlinear function of the 
radiative energy, i.e. p/CR = 1 x 10"£',r'-^. Last, we allow refine- 
ment with a criterion based on the radiative energy gradient. In 
each region where VEJE^ > 3 %, the grid is refined. 

Figure |3lfl) shows the radiative energy profiles at time t = 
1.4 X 10"^ for calculations run with a coarse grid of 16 cells and 
a maximum effective resolution of 512 cells (squares), and for 
calculations run with 2048 cells, taken to be the "exact" solu- 
tion (red curve). Because of the nonlinear opacity, the diffusion 
is more efficient in the high-energy region. The mean opacity at 
cell interface is computed using an arithmetic average, which 
is more adapted for the case of nonlinear opacity. The levels 
are finer (higher resolution) in high-radiative energy gradient re- 
gions. Note that we checked that we obtained similar results in a 
2D plane parallel case and in a 2D case with an initial step func- 
tion that maked an angle 7r/4 with the computational box axis. 
This validates our routine in the x and y directions. 

In Fig. [3/7) we show the evolution of the Li norm of the 
error as a function of the mesh spacing. The uniform grid points 
(diamonds) correspond to a calculations run with a number of 
cells ranging from 16 to 512 (i.e. ^ = 4 to ^ = 9) . When AMR is 
used (squares), the error is plotted as function of the minimum 
grid spacing, corresponding to effective resolutions ranging from 
32 to 512 cells. The coarse level remains unchanged, H^in - 4 
(i.e. 16 coarse cells). The advantage of the AMR is clear, the 
error remains identical compared to the uniform grid case, but 
the number of cells is greatly reduced (25 cells with ^'„,ax - 5, 38 
cells with /"max = 6, 59 cells with ^'^ax - 7, 83 cells with ^^ax - 8 
and 110 cells with ^'^ax = 9). The AMR implementation works 
well (second-order accuracy), and does not suffer from the fact 
that our scheme is only first order in space at the level interface. 
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which validates our scheme used to estimate the gradients at the 
cell interface in Sect. l3.7l Finally, as in the previous test, the error 
increases with h^, even if the diffusion problem is nonlinear for 
radiation. 



lution of the gas energy e, by solving the ordinary differential 
equation dTurner & Stonell200lh 




O.IOOOF 



0.0100 r 
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0.100 



de 

— = co-E, 
dt 



- 47TO-B(e). 



(22) 



We performed two tests, with two initial gas energies, e = 10'" 
erg cm"^ and e = 10^ erg cm"^. In both tests, the following 
quantities are taken constant: the radiative energy £1 = 1 x 10'^ 
erg cm"-', the opacity cr = 4 x 10"^ cm"', the density p = 10"^ 
g cm"-', the mean molecular weight /i = 0.6, and the adiabatic 
index y = 5/3. Figure |4] shows the evolution in time of the gas 
energy for the analytic solution (red line) and the numerical so- 
lution (squares). In the first calculations, where the initial gas 
temperature is greater than the radiative temperature, we used 
a variable time step At that increases with time, starting from 
10"^° s. This good sampling gives very good results. In the sec- 
ond case, we used a constant time step At = 10"'^ s. Although 
the sampling is bad at early times and longer than the cooling 
time, numerical solutions always match the analytic one. This 
validates our linearization of the emission term (arT'*) in equa- 
tion (fTsT i. 
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Fig. 4. Matter-radiation coupling test. The radiative energy is 
kept constant, Er - Ix I0'~ erg cm"^, whereas the initial gas en- 
ergies are out of thermal balance (e - 10^ erg cm"-' and e = 10"^ 
erg cm""*). Numerical (square) and analytic (red curve) evolu- 



Fig. 3. Nonlinear diffusion of an initial step function with AMR, 
the refinement criterion based on radiative energy gradients, (a) 
Radiative energy profiles at time t - 1 .4 x 10"^ (square - numer- tions of the gas energy are given as a function of time, 
ical solution, "exact" solution in red, run with 2048 cells). The 
AMR levels (right axis) are plotted in blue, (b) Li norm of the 
error as a function of li - 1/Ax, without AMR (diamond) and 

with AMR (squares), up to an effective resolution of 512 cells 4.4. w full RHD tests: Radiative shocks 
(the error is plotted as a function of the minimum mesh spac- 
ing, corresponding to the maximum resolution). The dotted line Testing the numerical method for radiative shock calcula- 
shows the evolution of the error as a function of h\ tions is a last important step that eveiy code attempting to in- 

tegrate RHD equations should perform (Haves & Norma nl2003l: 
Whitehouse & Bate 2006; Gonzalez et al. 2007). Following the 
Ensmanl ( ll994l) initial conditions, we tested our routine for sub- 
and super-critical radiative shocks. 

Initial conditions are as follows: uniform density po = 7.78 x 
10"'" g cm"-' and temperature To=10 K. The box length is 
L=7 X 10"'"cm, the opacity is constant (cr = 3.1 x 10"'" cm"'), 
ju = 1 and y = 7/5. The ID homogeneous medium moves with 
a uniform speed (piston speed) from right to left and the left 



4.3. Matter-radiation coupling test 

Another conventional test is the matter-radiation coupling. 
Consider a static, uniform, absorbing fluid initially out of ther- 
mal balance, in which the radiation energy E^ dominates and is 
constant. An analytic solution can be obtained for the time evo- 
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Fig. 5. Left: Temperature profiles for a subcritical shock with piston velocity v - 6 km s ' , at time t = 3.8 x lO'* s. Right: Temperature 
profiles for a supercritical shock with piston velocity v = 20 km s"', at time t = 7.5 x 10^ s. In both cases, the temperatures are 
displayed as a function of z - x - vt. The squares represent the gas temperature and the diamonds the radiative temperature. The 
red curves represent the gas and radiative temperatures obtained with a calculation using 2048 cells, which we take as the "exact" 
solution. The AMR levels (blue line - right axis) are overplotted. 



boundary is a wall. The shock is generated at this boundary and 
travels backwards. The piston velocity varies, producing sub- or 
super-critical radiative shocks. The AMR is used, and the re- 
finement criterion is based on the density and radiative energy 
gradients (30%), the grid has 32 coarse cells and we use five lev- 
els of refinement. We use the Minerbo flux limiter The time step 
is given by the hydrodynamics CFL for the explicit and implicit 
schemes. 

Figure |5] shows the gas and radiative temperatures for sub- 
and super-critical radiative shocks, as a function of z = x - vf, 
where v is the piston's velocity. The AMR is used in both cal- 
culations. The squares represent the gas temperature and the di- 
amonds the radiative temperature. The red curves represent the 
gas and radiative temperatures obtained with a calculation using 
2048 cells, which we take as the "exact" solution. The subcritical 
shock is obtained with a piston velocity v - 6 km s whereas 
the supercritical shock is obtained with v = 20 km s"'. In both 
tests, the occurrence of an extended, non-equilibrium radiative 
precursor is obvious. As expected, pre- and post-shock gas tem- 
peratures are equal in the supercritical case. 

For t he subcritical case, the postshock gas te mperature is 
given by (lEnsma nlll994tlMihalas & Mihalaslll984 



To 



2(7 -ly 

Kiy + 1)2 ' 



(23) 



where K - k/finiu is the perfect gas constant. For our initial 
setup, this analytic estimate gives T2 ~ 810 K. Numerical cal- 
culations give T2 ~ 825 K at time f = 3.8 x 10^ s, which agrees 
with the analytic estim ate comparable to val ues obtained with 
more accurate methods (iGonzalez et al.ll2007l) . The characteris- 
tic temperature T ~ 275 K immediately i n front of the shock 
agree s very well with the analytic estimate (iMihalas & MihalasI 
[1981 



276 K. 



(24) 



pvR V3 

This means that in front of the shock, the gas internal energy 
flux flowing downstream is equal to the radiative flux flowing 



upstream. The entire radiative energy is absorbed upstream and 
contributes to heat the upstream gas. Similarly, the spike temper- 
ature T+ ~ 1038 K also ag rees well with the analytic estimate of 
IMihalas & MihalasI (11984 



T2 + - — ^r_~ 980K. 
y + I 



(25) 



We note that the AMR scheme enables us to accurately 
describe the gas temperature spike at the shock. The medium 
around the spike is optically thin, and the numerical resolution 
in this region is therefore of crucial importance. The spike's 
amplitude varies according to the model used for radiation 
and to the effective numerical resolution. Thanks to the AMR 
scheme, the spike's amplitude is larger in the supercritical case, 
but not as large as thos e obtained with Ml o r VTEF models 
(iHaves & Normanll2003l iGonzalez et"ani2007h . However, this 
last test shows the ability of our time-splitting method to inte- 
grate the RHD equations. 

4.5. 3D dense-core-collapse calculations without rotation 

In this section, we perform calculations of a 1 Mq dense-core 
collapse without rotation, using our grey FLD solver We com- 
pare our FLD results for a model without in itial rotation with 
the ones obtained by iMasunaga et al.l ( 19981) and with ou r re- 
sults obtained with a ID code (see Commeryon et al.ll201 lb . We 
also qu a litativ ely c ompare our results with the p ioneered ones of 
iLarsonl (1 19691) and IWinkler & Newman! (Il980l) . This latter test 
provides a validation in 3D for star-formation purposes. 



4.5.1 . Initial conditions 

To facilitate the comparison with other studies, we used 
the same initial conditions as in ICommeryon et al.l (l2008l) and 
in Paper II and the Lax Friedrich Riemann solver We chose 
highly gravitationally unstable initial conditions. The initial 
sphere is isothermal. To - 10 K, and has a uniform density 
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Reference 






M 
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Tfc 


5,. 






(AU) 


(Mo) 


(Mo/yr) 


(C) 


(K) 


(K) 


(ergK-' g-') 




This work 


8 


2.1 X 10^^ 


3.7 X 10-' 


0.014 


396 


81 


2.11 X 10" 


24 


Masunasa et aL (1998^ 


~ 8 


~ 10-2 


~ 10-5 


0.002 


~ 200 


60 


2.08 X 10" 


6 


Commercon et al. (2011) 


7 


2.31 X IQ-^ 


3 X 10-5 


0.015 


419 


70 


2.02 X 10" 


19 


Larson 0969) 


^4 


~ 1 X 10-^ 






170 









Table 2. Summary of first core properties at time t =1.012 fff and = 2.7 x IQ-" g cm--'. Rf^, Mfc and Tf^ give the radius and 
mass of the first core, and the temperature at the first core border respectively. The mass accretion rate M and accretion luminosity 
^acc = GMfcM//?fc are also computed at the first core border and give the central te mperature and ent r opy, a ac.c is a typical 
accretion parameter The comparat ive values are roughly estimated at Pc ~ xlQ-" g cm-^ in lMasunaga et al.l (Il998h . at pc - IQ-"' 
gcm-^ in lCommer9on et all (1201 ih andatpc = 2 x 10-'° g cm-^ in lLarsonl(ll969h . 



po - 1.38 X 10-'^ g cm--'. The ratio a between initial thermal 
energy and gravitational energy is ~ 0.50. The initial radius is 
Ro = 7.07 X 10'* cm. The theoretical free-fall time is % = 57 
kyr. The initial isothermal sound speed is Cso ~ 0.19 km s-' and 
y = 5/3. The outer region of the sphere is at the same temper- 
ature as the core temperature, but is 100 times less dense. The 
sphere radius is equal to a quarter of the box length to minimize 

border effects. 

We use the set of opacities given by ISemenov et al. I (l200l 
for low temperature (< 1000 K), which we compute as a function 
of the gas temperature and density. For each cell we perform a 
bilinear interpolation on the mixed opacities table. Below 1500 
K the opacities are do minated by grain (silicate, iron, troilite, 
etc.). [Semenov et al.l (i2QQ3) take into account the dependence 
of the evaporation temperatures of ice, silicates, and iron on the 
gas density. We here use spherical composite aggregate particles 
for the grain structure and topology and a normal iron content in 
the silicates, Fe/(Fe + Mg)=0.3. 

4.5.2. Results 

To resolve the Jeans length, we use A^j = 10 (i.e. 10 
points per Jeans length) . iMasunaga et al.l (1 19981) showed that 
the first core properties are independent of the initial conditions 
for low-mass star form ation. We can then com pare our results 
with those obtained by IMasunaga et alj (Il998l) . even if we use 
different initial conditions. We also com pare our results w ith 
those obtained using a I D spherical code (lAudit et al.ll2002l) in 
ICommeryon et al.l (1201 ll) . 

Table |2] summarizes the first core properties obtained 
at time t =1.012 tg with RAMSES. The first core radius 
and mass are qualitatively similar to t he results obtained in 
other ID Lagrangean c alculations (see IMasunaga et al.l 119981 : 
ICommeryon et al.[ 1201 lb . even though we use a completely 
different hydrodynamical scheme (e.g. no artificial viscos- 
ity, Eulerian, etc.). The first core radiu s is ho w ever a fac- 
tor 2 greater than the o ne found in Larson I (119691) and 
IWinkler & NewmanI (Il980l) . who used simplified dust opacity 
models. Since the first core is mainly set by the opacity, this 
explains the differences. We define the first core radius as the ra- 
dius at which the infall velocity is maximal. The accretion rate 
on the first core is typical of low-mass star formation, ~ 10-^ 
Mo/yr We note that the value o-acc ~ 24 is relatively high, 
with Q-acc defined as M = cfaccC^Q/G (where Cso the isothermal 
sound speed). This indicates that our collapse mo del is closer to 
the dynamica l Larson-Penston collapse solution ( LarsonI 1969t 
|Penstonll969h than to the classical SIS model of lShuld 19771) . for 
which ttacc ~ 0.975. 

In Fig. |6| we show the profiles of density, radial velocity, 
temperature, optical depth, and integrated mass as a function 



of the radius and the temperature as a function of the density 
in the computational domain at time t -1.012 tg. All quantities 
are mean values in the equatorial plane. In the density profiles, 
all cells of the calculations have been displayed (blue points). 
The spread in the density distribution is very small. The spher- 
ical symmetry is thus well conserved in the 3D calculations 
w ith RAMSES. We compare these profiles with those obtained 
in lCommer9on et all (1201 lb . The density jump between the first 
core border and the centre is of the same order of magnitude 
as for the ID spherical case. The infall velocity at the shock is 
also comparable (~ 2 km s-'). The accretion shock takes place 
around t ~ 5 - 10, in the optically thick region. We do not see 
a jump in temperature through the accretion shock, which is a 
supercritical radiative shock. Eventually, we see from the tem- 
perature versus density plot that the thermal behaviour of the 
gas is not perfectly adiabatic in the central core. The first core 
is not fully adiabatic and is able to decrease its entropy level by 
radiating in the upstream material. The slight kink in the curve 
at r ~ 80 K (log(r) ~ 1 .7) corresponds to ice evaporation in the 
opacity table. The opacity decreases abruptly, this is the reason 
why the cooling is more efficient in that region. 



5. Summary and perspectives 

We have developed a full radiation-hydrodynamics solver 
using the flux-limited diffusion approximation, which is inte- 
grated in the AMR RAMSES code. Our solver uses a time-splitting 
integrator scheme and combines explicit and implicit methods. 
Each step of the integration is detailed in this work. The method 
was successfully tested in several conventional tests in ID and 
2D. We demonstrated that our method is second-order accu- 
rate in space, even when AMR is used. We also performed col- 
lapse calculations of a non-rotating dense core and successfully 
compa red our results with those o btained by ^M asunaga et al.l 
(Il998l) and lCommer9on et al.l (1201 ll) . which are based on differ- 
ent methods in ID spherical codes. Our method has thus been 
demonstrated to be robust and well suited for star formation. 
In Paper 11 we present detailed RHD calculations with a very 
high resolution of dense-core collapse in rotation. We showed 
that our method enables us to accurately handle the heating and 
cooling processes. Last but not least, we exte nded our method to 
the radiation-magnetohydrodynamics flows in lCommeryon et all 
(2010). 

The next step following this work will be to tune our solver 
for adaptive time-stepping to make use of all benefits of the 
AMR in RAMSES. For example, the next stages of the col- 
lapse, the second collapse and the second core formation, re- 
quire a huge amount of numerical resolution and the dynamical 
timescale becomes much shorter An adaptive time-step scheme 
is then suitable. Another improvement is to use a multi-group 
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Fig. 6. Profiles of density, radial velocity, temperature, optical depth, and integrated mass as a function of the radius and the temper- 
ature as a function of density in the 3D computational domain. All values are computed at time f =1.012 ts. 



approach in the radiation so lver Some attempts have b een pre- 
sented in the literature (e.g. IShestakov & Offnerll2008l) . but the 
computational cost remains too high nowadays compared to the 
grey model. 
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Appendix A: The super-time stepping versus the 
conjugate gradient 

In this appendix we present the super-time stepping (STS) 
method. It is used to solve parabolic equation systems, like the 
conjugate gradient (CG) we used previously. We implement the 
STS scheme into RAMSES. We compare the CG and the STS 
methods for the particular case of the ID linear diffusion test 
presented in ^4.11 

A.1. The super-time stepping 

The STS is a very simple and effective way to speed up 
explicit time-stepping schemes for parabolic problems. The 
method has been recently rediscovered in 'Alexiade s et al.l 
<!l996|), but it remains relatively unknown in computational as - 
trophysics (iMignone et"ani2007t lO' Sullivan & Downesll2006l) . 
The STS frees the explicit scheme from stability restrictions on 
the time-step. It can be very powerful in some cases and is easy 
to implement compared to implicit methods that involve matrix 
inversions. 

The STS is designed for time dependent problem such 

as 

^(t)+AU(t)^0, (A.l) 
at 

where A is a square, symmetric positive definite matrix. 
Equation lA.ll is rewritten with the corresponding standard ex- 
plicit scheme 

= U" - AtAU", (A.l) 

The explicit scheme is subject to the restrictive stability condi- 
tion 

p{l-AtA)<l, (A.3) 

where p( ) denotes the spectral radius. The equivalent CFL con- 
dition is 

2 

At < Afexpi = , (A.4) 

^max 

where /Imax stands for the highest eigenvalue of A. For the ID 
heat equation du/dt - x^^^^ discretized by standard second- 
order differences on a uniform mesh, we have Amax - 
(Afexpi = Ax2/2a-). 

In the STS method, the restrictive stability condition is 
relaxed by requiring the stability at the end of a cycle of A^sts 
time-steps instead of requiring stability at the end of each time 
step Af. It leads to a Runge - Kutta- like method with A^sts stages. 
Following lAlexiades et al.l (Il996h . we introduce a superstep 
AT - Yj^^i consisting of A^sts substeps ti,T2,- ■ •,TNa^- The 
idea is to ensure stability over the superstep Ar, while trying to 
maximize its duration. The inner values, estimated after each 
Tj, should only be considered as intermediate calculations. Only 
the values at the end of the superstep approximate the solution 
of the problem. 
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The new algorithm can be written as 

/'A'sts 



Y]{l-TjA) 



U", 



and the corresponding stability condition is 

< 1. 



Y]{l-TjA) 



V./=i 



(A.5) 



(A.6) 



In order to find AT as high as possible, the properties of 
Chebyshev polynomials are exploited, providing a set of opti- 
mal values for the substeps given by 

VstsjCOS 



Tj - Afexpl 



(-1 



(A.7) 



where Vjts is a damping factor that should satisfy < Vjts < 
'tmin/'tmax- The supcrstcp AT" is given by 



Ar 



"SIS 



At 



expl 



o 1/2 



(1 



sts ^ 



(i+y:52^+(i-v:52^» 



(A.8) 

Note that Ar — > A^^(j,Afexpi as Vsts 0. The method is unstable 
in the limit Vjts = 0. The STS method is thus almost A^sts times 
faster than the standard explicit scheme. When Ar is taken to 
be the advective (CFL) time step At while coupling with the hy- 
drodynamics, the STS requires only approximately ( Af/Afexpi)'^^ 
iterations rather than Af/Afexpi with an explicit scheme. 

A.2. The STS implementation for fhe FLD equation 

The STS scheme replaces the implicit radiative scheme pre- 
sented in Sect. 13.41 Equations of system (fT4l i written with an 
explicit scheme become 



C,T"+' - C,T" 



At 



-E" 



cA" 



• 4p"c(«R(n'-£") 



At K'^p" 

(A.9) 

The explicit time step Afgxpi is estimated using values at time 
n. The next step consists of estimating values of A^sts and Vsts, 
the latter dep ending on the spectral p roperties of A. However, as 
mentioned in lAl exiades et alJ (Il996h . it is not required to have 
a precise knowledge of the spectral properties for the method 
to be robust. A^sts and Vsts are thus arbitrary chosen by the user 
Instead of executing one time step of length Afgxpi, one executes 
supersteps of length AT". A^sts substeps ti,T2, ■ ■ ■,TNsts are thus 
performed without outputing until the end of each superstep. 
When the STS is coupled to the hydrodynamics solver, the cy- 
cle is repeated until the time step, given by the hydrodynamical 
CFL condition, is reached. Superstep Ar and substeps t, are re- 
estimated at the end of each cycle. 

A.3. Comparison w/fh fhe conjugate gradient metliod 

To compare the STS with the CG algorithm we used through- 
out, we consider the test case presented in Sect. 14.11 The equa- 
tion to integrate is simply 

dt \ 3pKji^ 



-V£, 



0. 



(A. 10) 



The initial setup is identical to those in Sect. 14.11 It consists of 
an initial pulse of radiative energy in the middle of the box. We 
present here calculations made with either the STS method or 
the CG algorithm. In both cases, CG and STS are applied over 
an arbitrary time step Af that simulates the time step that would 
be given by the hydro CFL. All calculations were performed on 
a grid made of 1024 cells. In the STS calculations, for each value 
of Af, calculations have been performed using various values of 
A^sts and Vsts- For the CG method, only the convergence criterion 
econv chan ges. 

Figure lA.ll shows the radiative energy profiles at time f = 
1 X 10 s for all calculations we performed. In all panels, the 
analytic solution is plotted (black line). The two upper plots give 
results for the CG method. For Af > lO"''*, the accuracy is very 
limited. We also see that for Af > 10"'^, the diffusion wave does 
not propagate at the right speed. The total energy is conserved, 
but the diffusion wave does not have the correct extent. But the 
STS results are much more accurate, except for A^sts - 20 and 
Vsts = 1 X 10"^. By construction, STS is expected to be more 
accurate. The stability is poor when A^sts - 20 and V sts = 1 X 10~^ 
because Vsts is close to the stability limit (see lAlexiades et alJ 
Il996h . 




1.0 



Fig. A.2. Comparison of calculations done using STS or CG and 
a variable time step given by Af = 1 x 10"'^ * 1.05"*"^p, where 
istep is the index of the number of global (hydro) time steps. 
Results are given at time f = 1 x lO^'-'s. 



In Table lA.ll we give the CPU time and A^itei , which corre- 
sponds to the number of iterations for the CG and to the number 
of substeps for the STS. The number of operations per iteration 
in the CG and per substep in the STS are equivalent, since it 
involves the same number of cells (1024). The CPU time spent 
with the STS is ten times smaller than the one of the CG method. 
The STS also requires often twice less iterations than the CG. 
The bottom lines give the results for calculations made with a 
variable time step, which increases with time. The corresponding 
profiles are plotted in Fig. IA.2I The STS remains more accurate 
in this case than the CG, which is quite accurate over more than 
three orders of magnitude. The CG gives good results, because, 
thanks to the variable time steps, the diffusion wave propagates 
at a correct speed. Indeed, at f = 0, the gradient of radiative en- 
ergy is steep and the diffusion wave speed is very high. Using an 
initial short time step Af = 10"'^ enables us to be closer to the 
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Method 


Parameters 


At 


CPU time (s) 


A'itcr 


CG 


F - 1 X 10"'' 

*=t;onv — i ^ iu 


1 X 10""' 


82.9 


10805 


CG 




1 X 10"'* 


68.5 


14623 


STS 


V - 001 N - 10 


1 X 10"" 


20.4 


9000 




, — 1 X 10"'' 


1 X 10"'^ 


21.2 


47^8 

^ 1 JO 




, — 1 X If)"** 


1 X 10"'^ 


27.8 




STS 


= 1 X 10"'' = 20 


1 X 10"'^ 


2.8 


2107 


STS 


V . - 001 N . - 5 


1 X 10"" 


2.9 


2408 


STS 


„ _ f) nni N . - 90 


1 X 10"" 


2.9 


2408 


STS 


V . - 001 N . - 10 


1 X 10"" 


2.9 


2408 


STS 


V . - 005 N . - 5 


1 X 10"" 


3.1 


2709 


STS 


V . - 005 N . - 1 


1 X 10"" 


3.04 


2709 




, — 1 V If)"'' 


1 X lO"'-* 

1 A IL/ 


1 1 A 


9848 


cc. 


, — 1 V If)"** 
tconv — 1 A lU 


1 X 10"'* 
1 A lU 


1 S A 


^8Q9 
joy A 


STS 


V , = 1 X 10"'' N , = 20 


1 X 10"'* 


0.6 


600 


STS 


V . - 001 N . - 5 


1 X lO"''* 


0.99 


1470 


STS 


V . - 001 N . - 70 


1 X lO"'"* 


0.75 


900 


STS 


„ _ f) nf)| N - 10 


1 X 10"'"' 


0.69 


780 


STS 


V . - 005 N . - 5 


1 X lO"'"* 


1.1 


1620 


STS 

O 1 o 


, , _ n fin's N —in 


1 X 10"'"* 


87 


1 170 


cc 


, — 1 V If)"'' 
tconv — 1 A lU 


S X 1 0"'* 

J A LV 


v.y 


1 


CC 


f: — 1 V If)"** 
tconv — 1 A lU 


J A L\J 




z.j\jj 


STS 


V , = 1 X 10"'' N , = 20 
'^StS 1 /\ , lists 


5 X lO"''' 


0.39 


390 


STS 


V - n nni n — 

^^sts — U.UUi , lists — J 


5 X lO"'"* 


8 


1326 


STS 




5 X 10"'** 


0.56 


756 


STS 


V . - 001 N . - 10 

•'^sts — U.UUl , lists — 


5 X 10"'"* 


0.46 


534 


STS 


V - OOS N — S 

t^sts — U.UUJ , lists — 


S V lO"'"* 

J A 1\J 


87 


1476 


STS 


„ _ n nn^ ]\j -10 

^sts — vJ.UUJ , iNsts — 


5 X 10"''' 


0.68 


1032 


CC, 


, — 1 V 10"'' 

tconv — 1 A lU 


1 X 10"'' 

1 A lU 


A f, 


1 1 ^S 

i LJJ 


cc 


t _ 1 y 1 n-8 

tconv — 1 A lU 


1 V 10"'-^ 
1 A lU 


J.O 


1 ^QQ 
Ljyy 


STS 


Vsts = 1 X 10"'' , Nsts = 20 


1 X lO"''' 


0.36 


351 


STS 


Vsts = 0.001 , Nsts = 5 


1 X 10"'3 


0.77 


1311 


STS 


Vsis = 0.001 , Ns,s = 20 


1 X 10"" 


0.52 


729 


STS 


Vsis = 0.001 , Ns,s = 10 


1 X 10"'^ 


0.43 


495 


STS 


Vsts = 0.005 , Nsts = 5 


1 X 10"'3 


0.83 


1467 


STS 


Vsis = 0.005 , Ns,s = 10 


1 X 10"'^^ 


0.65 


1020 


CG 


tconv = 1 X 10"** 


1 X 10""= * 1.05""=!^ 


19 


4680 


STS 


Vsis = 0.005 , Ns,s = 10 


1 X 10"'* * 1.05'"'='^ 


1.7 


1782 



Table A.l. Summary of calculations plotted in Fig. lA.ll CPU time, the number of iterations in the CG method or the number of 
substeps in the STS method are given for various time steps and various values of tcotiv for the CG, and Nsts and Vsts for STS. 



CFL condition associated to the diffusion wave speed. Then, ra- 
diative energy gradients and the former CFL condition relax and 
the time step can increase with time. This relaxation on the inte- 
gration time step enables us to maximise the accuracy of implicit 
methods using a subcycling scheme based on the diffusion wave 
speed propagation. However, this speed remains quite difficult to 
estimate. 

Eventually, we must conclude by pointing out that even if 
the STS method is well adapted for this problem, it remains very 
limited for star-formation calculations. Indeed, the diffusion time 
is very short compared to the dynamical time estimated as the 
free-fall time (see Fig. lA.3b and then, the STS requires too many 
substeps. The convergence of the CG depends on the nature of 
the problem and is not affected by strong differences between 
the diffusion and the dynamical times. Moreover, we never en- 
counter these steep gradients in the radiative energy distribution 
in star-formation calculations. The STS could be efficient only 
within the fragments, where the diffusion time is very long. This 
is the reason why we only use the CG method throughout. An 
alternative but non-trivial solution would be to couple the CG 
and the STS methods. 
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Fig.A.3. Contours in the equatorial plane of the ratio between 
diffusion and free-fall times for collapse calculations. The diffu- 
sion time is estimated as Tdiff = ^r^, where / is the local Jeans 
length. 
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Fig. A.l. Comparison of the numerical solutions using STS or CG with the analytic one (black line) at time f = 1 x 10 s. The 
color curves depict numerical solutions obtained with timestep Af equals to 1 x lO"'-' (blue), 5 x 10"'"^ (red), 1 x 10"''' (green), 
1 X 10"'^ (yellow) and 1 x 10""' (cyan). 



